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PREFACE 


In  this  Memorandum  the  authors  present  further 
mathematical  results  from  their  study  of  radiative 
transfer,  using  invariant— imbedding  techniques.  This 
subject  has  important  implications  for  meteorology, 
astrophysics,  and  nuclear-blast  detection. 

The  research  presented  here  was  sponsored  by  the 
Advanced  Research  Projects  Agency. 


SUMMARY 


In  the  study  of  radiative  transfer,  the  classical 
Boltzmann— type  approach  focuses  attention  on  the 
statistical  behavior  of  the  particles  within  the  medium. 
There  are,  however,  analytical  and  computational 
advantages  in  the  invariant— imbedding  approach,  which 
leads  to  the  determination  of  the  diffusely  reflected 
and  transmitted  intensities  in  slab  geometry.  In  the 
present  paper  it  is  shown  that  the  invariant— imbedding 
technique  is  also  extremely  useful  for  the  study  of  the 
Internal  intensity  of  radiation  in  a  finite,  homogeneous, 
flat  layer.  In  other  words,  the  technique  provides  us 
directly  with  a  complete  solution  of  scattering 
problems  of  light  in  a  finite  plane-parallel  atmosphere. 
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INVARIANT  BIBEDDING  AND  A  RESOLVENT  OF 
THE  PHOTON-DIFFUSION  EQUATION 

1,  INTRODUCTION 

While  it  is  known  that  the  invariant-imbedding 
technique  is  powerful  for  the  study  of  the  diffuse 
reflection  and  transmission  problems  of  radiation  in 
slab  geometry,  the  usefulness  for  inquiring  into  the 
outward  and  inward  intensities  of  radiation  at  any 
level  in  the  medium  is  not  so  well  -understood  (see 
Bellman  and  Kalaba  [  1J  ;  Wing  [  123  ). 

Recently,  Sobolev  ([?3,  [83  )  has  shown  that  when 
the  Milne  type  of  a  linear  integral  equation  with  an 
arbitrary  source  distribution  is  considered  in  the  case 
of  isotropic  scattering,  the  solution  is  obtained  by 
reducing  the  determination  of  the  resolvent  in  terms  of 
two  optical  depth  arguments  of  the  integral  equation  to 
that  of  the  function  $(t)  depending  only  on  a  single 
argument  r.  The  partial  differential  equation  of  the 
resolvent  was  found  by  using  the  physical  and  probabilis¬ 
tic  manipulation.  The  invariant-imbedding  technique  was 
applied  to  the  solution  of  an  integral  equation  for  the 
emission  probability  distribution  in  a  semi— infinite, 
homogeneous,  flat  layer,  and  then  the  resolvent  of  the 
photon-diffusion  equation  was  expressed  in  terms  of  the 
fundamental  function  ®(t)  in  connection  with  the 
integro-differentlal  equation  of  the  emergence 
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probability  of  a  photon  (see  Ueno  [93).  Furthermore, 
the  invariant— imbedding  approach  was  used  to  determine 
the  radiation  field  within  the  high— temperature  stagnation 
region  ahead  of  the  blunt-nosed  model  in  a  hypersonic, 
rarefied  gas  stream. 

In  this  paper  we  shall  show  how  the  invariant- 
imbedding  technique  can  be  used  to  find  the  source 
function,  in  the  radiation  field  with  arbitrary  distri¬ 
bution  of  emitting  sources,  by  the  agency  of  the 
resolvent  of  the  photon— diffusion  equation.  Thus,  by 
evaluating  the  requisite  intensity  inside  the  homogeneous 
medium  according  to  the  standard  procedure,  a  complete 
solution  of  the  transfer  equation  of  radiation  is  obtained. 
Finally,  the  relationship  between  the  scattering  function 
and  the  resolvent  is  given  by  means  of  the  emergence 
probability  of  photon. 

It  is  of  interest  to  mention  that  the  principles  of 
invariance  can  also  be  applied  to  find  the  intensity  at 
any  level  t  as  follows  (see  Chandrasekhar  [43;  Horak  t 53 ) 


i(t,-hjl) 


-  +gr. 


s  (  t2— ■ r;  [I , \1 » )  I  (  t ,-+1 » )  dp, f 


0 


It(t,-+Jl)  +  ^  J  S(T;p,p«)l(^^,)ddS 


0 


where  IT  _T(T^+M')  (0  <  \x  <  1)  is  the  intensity  which 

would  be  emitted  in  the  outward  direction  at  any  level  if 

* 

there  were  no  layer  between  0  and  t,  and  I  (t,-jj.) 

ft 
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(0  <  |x  <  1)  is  that  in  the  inward  direction  at  any  level 
if  there  were  no  layer  from  t  to  Provided  that 

I*^_t(t,-hx),  I*(T,-JX),  S(t;ix,ix0),  and 

are  known,  the  procedure  of  successive  approximation 
leads  to  the  requisite  internal  intensities.  However, 
the  approximate  values  of  I(t,hhx)  for  starting  the 
procedure  should  be  evaluated  by  solving  approximately 
the  transfer  equation.  In  the  one-dimensional  case  it  is 
shown  that  the  reflected  and  transmitted  fluxes  can  be 


considered  as  fundamental  functions  from  which  all  other 
functions  can  be  derived  (cf.  Bellman,  Kalaba,  and  Wing 
[23 ).  Then  the  present  approach  may  be  called  a  direct 
procedure  for  obtaining  the  internal  flux,  whereas  the 
above  procedure  is  an  indirect  one. 

The  powerfulness  of  the  invariant— imbedding  technique 
for  obtaining  a  complete  solution  of  the  transfer  problems 
of  radiation  in  a  finite  flat  layer  lies  in  the  fact  that 
the  back— scattering  process  does  not  play  a  role  in  the 
determination  of  the  radiation  field  in  slab  geometry 
(see  Uesugi  [ 111  ) . 

Later,  the  present  approach  may  be  applied  to  the 
determination  of  the  perturbation  velocity  distribution 
in  the  linearized  version  of  the  Boltzmann  transport 
equation,  based  on  the  Krook  kinetic  model  in  the  classical 
Couette  flow  problem. 


2.  MILNE1 S  FIRST  INTEGRAL  EQUATION 


Consider  the  diffuse  reflection  and  transmission  of 
parallel  rays  by  a  homogeneous  flat  layer  of  finite 
optical  thickness  —  t  >  0,  scattering  radiation 
isotropically. 

Let  a  parallel  beam  of  radiation,  of  net  flux  ttF 
per  unit  area  normal  to  itself,  be  falling  on  the  surface 
t  =  tq  of  the  medium  in  the  direction  —  M-0  (0  <  M-q  <  1), 

where  m-q  denotes  the  cosine  of  the  inclination  to  the 
inward  normal  at  t  =  tq.  Following  the  notation  of 
Chandrasekhar  [  4] ,  the  equation  of  transfer  appropriate 
to  the  present  case  is  written  in  the  form 


(2.1) 


,  dl(T,p) 

dr 


I(t,h) 


I(t,U'  )dp* 


Fe 


-('T-Tj/d 


3 


where  is  the  albedo  for  single  scattering,  i.e.,  the 
probability  of  survival  of  a  photon. 

Consider  a  photon— diffusion  process  in  terms  of  the 
emission  probability  p(|a,t;t0,t1)  that  a  photon  absorbed 
at  the  level  t  (tq  <  T  <  will  reappear  in  the 

direction  +  \i  (0  <  p.  £  1)  in  the  radiation  emerging 
from  the  surface  t  =  For  convenience,  we  suppress 

the  parameters  and  in  p(ii,,T;T0,r1) . 

Taking  into  account  the  random  character  of  the 
photon— diffusion  process,  an  integral  equation  governing 
p(m,,t)  is  expressed  in  the  form  (see  [9]) 


(2.2)  [l  -  cd/l!  T{p(p,t)}  =  cue 


-(t— Tn)/jX 


where  1  on  the  left-hand  side  is  an  identity  operator, 
and  the  truncated  Hopf  operator  T  is  given  by 


(2.3)  AT(f(t)}  =  |  f  Tl  f(t)E1(|t  -  t| )dt. 


In  Eq.  (2.3),  is  the  first  exponential  integral  for 


the  positive  real  argument 


(2.4) 


E1('r)  = 


-t/U  dp. 


With  the  aid  of  the  combined  operations  method  (see 
Busbridge  L3^)*  it  is  readily  shown  that  Eq.  (2.2)  is 
equal  to  the  auxiliary  equation,  i.e.,  the  Milne  first 
integral  equation  for  the  source  function  ^(^p). 

Furthermore,  consider  the  scattering  of  light  in  an 
emitting  atmosphere  of  finite  optical  thickness  —  tq 
whose  distribution  of  emitting  sources  is  given  by  B(t), 
The  equation  of  transfer  appropriate  to  this  case  takes 
the  form 


(2.5) 


=  i*(t,p)  -^(t). 


where  the  source  function  5{t)  is 


r(T)  =  §  f  i^t^p1)^1  +  b(t) 


(2.6) 


* .  . 

Let  I  (T,4p),  0  <  p  <£  1,  be  the  intensity  of 

radiation  at  level  t  directed  towards  the  surface 
t  =  and  let  I*(t,~p),  0  <  p  <  1,  be  that 

directed  towards  the  surface  t  =  t^. 

Then,  as  the  formal  solution  of  Eq.  (2.5)  we  have 

(2.7)  I*(t,-hi)  =  J*1  5(t)e"(t^^ 

T 

(2.8)  =  /  T 

To 

Equations  (2.7)  and  (2,8)  show  that  once  ^(t) 
has  been  found,  the  inward  and  outward  intensities  at 
each  level  can  be  determined. 

In  what  follows,  it  is  shown  that  a  resolvent  of 
the  photon— diffusion  equation  permits  us  to  evaluate  the 
source  function  S'(t). 

The  Milne  first  integral  equation  for  ^(t)  is 
given  ( see  t  3^  )  by 

(2.9)  tl  -  oiKlT{s'(t)}  =  B(t)  . 

With  the  aid  of  the  resolvent  K(t,  t;  t0,ti  )  *  E0.» 

(2.9)  takes  the  form 

(2.10)  y(r)  =  B(t)  +  J  1  K(T,tjTQ,T1)B(t)dt, 

T0 

where  the  resolvent  K(T,tjTQ,T;L)dTdt  is  the  probability 
that  a  photon  emitted  in  the  optical  depth  interval 
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(t,t  +  dt)  will  be  re-emitted  between  the  optical  depths 
t  and  t  +  dT  after  one  or  more  scattering  processes. 
Then,  given  the  resolvent  K— function,  Eq,  (2.10)  enables 
us  to  evaluate  the  source  function  (see  Fig.  l). 


Incident  intensity 


t  —A  T 
0^  0 


Reflected  intensity 


Transmitted  intensity 
T  t  ^ 


Fig.  1  Reflection  and  transmission  of 
radiation  in  rod  geometry 


3.  RESOLVENT  OF  THE  PHOTOI'L-DIFFUSION  EQUATION 

Let  0(x;Tq,t^)  be  the  probability  that  a  photon 
absorbed  at  t  =  x  (T0  <  x  <L  Ti[)  will  reappear  in  the 
radiation  emerging  from  the  surface  t  =  after 
successive  scattering  processes. 

On  replacing  t  in  Eq.  (2.2)  by  x,  multiplying 
it  by  l/(2[x)  and  integrating  with  respect  to  p  over 
(0,1),  we  have 

(3.1)  §(x;TQtT1)  =  g(x-TQ)  +J1  cE)(xt;T0,T1)k(x,,x)dx' , 

T0 

where  tq  <  x  < 


0 


g  (x  Tq  )  o  ^  E1(x-T0),  k(x»,x)  =  §  E2  (  I  X 1  —X  1  )  . 


(3.2) 
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Altematively, 


(3.3)  ^(xiTQ,^)  =  g(x-Tn)  +  /  K^x';^,-!:.,  )g(x'-Tn)dx» , 


‘O'  1J 


where  K(x,x* is  the  resolvent. 

It  is  manifest  that  in  Eqs.  (3.1)  and  (3.3)# 
k(x,x* )  and  K(x,x’;Tq,t^)  together  are  symmetric  with 
respect  to  x  and  x*  (see  Mikhlin  [6]). 

Adding  a  layer  of  infinitesimal  optical  thickness 
A  to  the  surface  t  =  tq  of  the  medium,  we  now  consider 
the  invariant— imbedding  procedure. 


(3.4)  qp(x)  =  g(x-Tn+A)  +  f  1  cp(x*  )k(x*  ,x)dx* . 


rQ-A 

On  the  other  hand,  it  is  well  known  that  the 
resolvent  fulfills  an  integral  equation  given  [ 6]  by 


(3.5)  K(x,y;TQ,T1)  =  k(x,y)  +  f  1  K(x,y‘ ; Tn,Tn )k(y» ,y)dy* . 


'O'  1' 


Then,  making  use  of  Eq.  (3*5)#  Eq.  (3.4)  reduces  to 


(3.6)  <p(x)  =  g(x+A-Tn)  +  9(t_— 6A)k(Tn— 0A,x)A 


K(x,X* i T0,T1)g(x*+A— Tn)dx’ 


+  cp('T0-e0A)A  f  1  K(x,y;Tn,T1  )k(Tn-enA,y)dy 


0'  ‘1'^  'CT'0L 


+  6( A), 


where  0  <  0  (or  6^)  <  1.  The  details  of  the  compu¬ 

tation  will  he  referred  to  in  the  Appendix. 

On  the  other  hand,  recalling  Eq.  (3.3) ,  we  get 

(3-7)  <P(x)  a  g(x-TQ+A)  +J1  K(x,x«  jT0-A,T1)g(x,-T0+A)dx' 

“  g ( x+A  tq )  +J1  K(x,x' ^^^(x'-TQ^dx1 
T0 

-A  J  1  (x,x*;T0,T1)g(x!-T0+A)dx» 

To  ° 

+  K(x,T0-eAJ*T0-A,T1)g(A(l-e))A  +  0(A), 


where  0  <  0  <  1. 

With  Eq.  (3.3)*  the  combination  of  Eqs.  (3.6)  and 

(3.7)  provides  us  with 

(3.8)  9(r0-eA)k(T0-eT,x) 

+  cp(T0-e0A)  J  1  K^x^T^T^k^-eQ^x'  )dx» 

To 

=  K( x, tq  0A j  Tq  A, r x ) 

*  cp(To-0A)  -J  1  K(tq— 6A,xf ;  tq— A,T1)g(x‘— TQ+A)dxf 


1  dK 


(x,X'iT0^i 


'0 


Passing  to  the  limit  as 


)g(x‘— TQ+A)dx'  +  0(A). 
A  0,  from  Eq.  (3.8) 


we 


have 


Furthermore,  in  what  follows,  we  shall  derive 
another  type  of  partial  differential  equation.  In  a 
manner  similar  to  that  given  before,  we  shall  add  a 
layer  of  small  optical  thickness  A  to  the  surface 
t  «  t1. 

Let 


(3.10)  $(x)  =  g ( x— t 0 )  +  /  1  f(yj  )k(x*  ,x)dx* . 


T. +A 


'0 


Then,  allowing  for  Eq.  (3.5),  Eq.  (3.10)  reduces  to 


(3.11)  Tjr(x)  =  g(x-TQ)  +  A)k(x1+  A,X)A 

+  f  1  K(x,x,;T0,T1)g(x,-T0)dx*  +  ^(T1+e0A)A 
T0 

*  J  1  K(x,x* ;T0,T1)k(T1+e0A,x* )dx  +  0(A), 

T0 

where  0  <  6  (or  6Q)  <  1. 

On  the  other  hand,  we  have 

(3.12)  f(x)  =  g(x-TQ)  +J1  K(x,x,;T0,T1)g(x'-T0)dx‘ 

T0 

P 

+  a  J  -&T  (x,y;T0,T1)g(y-T0)dy 

To  1 

+  K(x,T1+eA;T0,T1+A)g(T1+eAr-T0)A  +  0(A), 


where  0  <  6  <  1. 
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Making  use  of  Eq.  (3*3)  and  combining  Eqs.  (3.11) 
and  (3.12),  we  get 

(3.13)  ^(T1+eA)k(T1+0A,x) 

+  ^(T1+eQA)A  J  1  K(x,xt;T0,T]L)k(T1+e0A  jX'Jdx* 

T0 

=  K(x,t1+0A;to,t1+A) 

/■  jr  r|-i^^ 

1  K(T1+eA,y;T0,T1+A)g(y-T0)dy 

L  To 

+  f  1  (^.y;'r0^1)g(y-T0)dy  +  ^(a). 

To  1 

Using  Eq.  (3.1)  and  putting  A  ->  0,  from  Eq.  (3.13) 
we  have 

(3. 14)  (x,y;  Tq,^)  =  K^T^T^T^K^y;  T^)  . 

Taking  into  account  the  optical  homogeneity  of  the 
medium  and  putting  tq  =  0,  Eq.  (3-14)  becomes 

(3.15)  (x,yjr1)  =  K(T1-x,0jT1)K(0,T1-yjT1), 

which  is  equivalent  to  that  yielded  by  Sobolev  [83. 
Because  of  the  homogeneous  medium,  we  get 

(3.16)  K(x+A,y+AjTQ,T1)  =  K(x,yiT0~A,T1-A) . 
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Then, 


,,  17.  c*K  .  BK  .  3k  .  BK 

(3.1?)  BIq  +  B5E  +  By  *  Bx^ 


=  0. 


On  recalling  Eqs.  (3*9)  and  (3.1^)#  from  Eq.  (3*17) 
we  obtain 


(3.18)  |f  +  |p  -  K(x,T0;T0,T1)K(T0,y;T0,T1) 

—  K(x, J Tq, )K( ,y j Tq, ) . 

When  we  allow  for  the  homogeneity  of  the  medium  with 
respect  to  the  optical  depth,  and  write 

(3«19)  K(x,TqJ  Tq,T^  )  =*  K(x— TQ,OjTj — Tq), 

K(x,T^j  Tq,T^)  =  ji ( T^—X,  Oj  Tq)  , 

Eq.  (3.18)  becomes 

(3.20)  |§  +  ||  -  K(x^r0,0;T1-t0)K(o,y-^0;T1-^0) 

-  K(  t2-x,  Oj  'fj— tq)K(  0,  Tx-yj  ) . 

Thus,  for  x  >  y  we  get 


(3.21)  K(x,y;  Tq,^)  =  KCx-y^j^-TQ) 


|K ( t+x— y— tq , 0 ;  t1~tq )K( 0, t— tq; t^-Tq ) 
K(  Tf-t-x+y,  0;  t1-t())K(  0,^-t;  ^-Tq) 


dt. 
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Then,  the  function  K(x,y;TQ,T1)  of  the  two 
arguments,  x  and  y,  Is  expressed  in  terms  of  the 
function  K(x— of  a  single  argument  x. 

Mien  t0  *=  0,  Eqs.  (3-20)  and  (3*21)  reduce  to  those 
given  by  Sobolev  [83,  with  the  aid  of  the  probabilistic 
approach. 

From  Eq.  (3*5)  we  have 

(3.22)  KU-TqjOjt^q) 

=  g(x-T0)  +J1  K(t-T0,0;'r1-T0)k(t,x)dt. 

T0 

Thus,  given  g(x— ■ tq),  the  one-parametric  resolvent 
K(x— Tq, 0; Tj—Tq)  can  be  determined  from  Eq.  (3.22). 

In  a  manner  similar  to  that  given  by  the  authors  in 
a  preceding  paper  [93,  we  shall  verify  an  identity 
K(x— Tq,0; Tj—Tq)  =  *(x;Tq,t^).  Following  the  notation  of 
Chandrasekhar  [43,  X(p;T0,T1)  is  given  by 

(3*23)  X()J.jT0,T1)  =  1  +  jy  J  S(Tq,T^ Jli.,11  * ) 

o'" 

where  the  scattering  function  S  is 

r  Ti  ~(t~T  )/\i 

(3.24)  S(t0,t1;i.l,m.0)  »  J  p(p0;t)e  u  dt. 

T0 

Putting  t  =  tq  in  Eq.  (2.2)  and  recalling  Eq. 
(2.10),  we  have 

(3.25)  P(y-,T0)  =  05  +  CD  J  1  K(TQ,tjT0,T1)e  0  dt. 
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Writing 


(3.26)  p(h,TQ)  =  oSC(ix;t0,t1), 
from  Eq.  (2.25)  we  get 

(3.27)  XCixjTq,^)  =  1  +  JTl  T°  K(T0,t+T0;T0,T1)e~t>^dt. 

0 


On  the  other  hand,  Eq.  (3.23)  becomes 

(3.28)  X(H;T0,T1)  =  1  +  J  1  °  ^(t+TQjT^T^e^/^dt, 

0 

Recalling  Eq.  (3.19),  the  comparison  of  Eqs.  (3.27) 
and  (3.28)  yields 

(3.29)  K(x-T0,OjT1-TQ)  a  ^(xjTqjT-jJ. 

Then,  Eq.  (3.21)  for  x  >  y  becomes 

(3-30)  K(x,y;T(),T1)  =  ©(x-y;^,^) 


'0 


® ( t+x-y-x0 j tq, r2 ) 0 ( t— tq j tq , r1 ) 


[  Tx-t-x+y  j  t0,t1)*( Tx-t ;  tq , t 2 ) 


Recently,  expanding  the  function  <I>(x;0,t1)  by  a 
power  series  of  co 

(3.31)  <t>(x;0,T.)  =  2  to1*,  (x,t .), 

x  i=l  1  1 


dt. 


Sobolev  and  Minin  [8]  computed  numerically  3>i(x,T1)  for 
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several  values  of  They  mention  that  the  approxima¬ 

tion  of  the  above  result  for  the  computation  of  the 
resolvent  is  of  practical  interest  only  when  just  a  few 
terms  of  the  expansion  formula  are  sufficient  for  the 
accuracy  required. 


APPENDIX 


Allowing  for  Eq.  (3.5),  recalling  the  symmetrical 
property  of  K(x,x' ;  t0*ti)  k{x,x*)  with  respect 

to  x  and  x1,  and  inverting  the  order  of  integration, 
we  have 


(A.l) 


'0 


1  K(x,x‘;'r0,T1)dx'  J  1  9{y)k(y,x‘  )dy 

T0 

ri  r  Ti 

M  k(x,x,)<3xl  /  9(y)k(y,x’  )dy 


0 

+  /  1  dx»  J  1  K(x,y;T0,T1)k(y,xI  )dy 
T0  T0 


fpCyOkCySx1  )dy 1 , 


'0 


and 
(A. 2) 


SCx'+A—Tq)  +  cp(Tn—  0nA)k(Tn—  enA,X*  )A 


Q  *Q  '0<- 


dx* 


'0 


'0 


'0 


0 


K(x,y,*T0,T1)k(y,x«  )dy 


cp(x' )  —  J  cp(x"  )k(x"  ,x'  )dx"j  dx1 
T0 


K(x,y; T0,T1)k(y,xt )dy 


(x,yj Tn,r 


O'  i 


cp(x*  )k(x*  ,y)dx* 


'0 


'0 
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T  o  T  P  T 

x6x'  x  cp(x"  )k(x,!  ,x'  )dx"  '  a 


K(x,y; T0,T1)k(y,x* )dy 


T 

1 


0  '0  ’0 
T-, 


k(x,x,)dxt  J  9(y)k(yjX* )dy. 

To  To 

Making  use  of  Eqs.  (A.l)  and  (A. 2),  Eq.  (3.4)  becomes 


(A. 3)  <p(x)  =  g  ( x-f-A—T  q  )  +  cp  (  t  n~6A )  k  ( rn-6A,  x )  A 


0 


'0 


cp(x’  )k(x’  ,x)dx'  +  {9(a) 


'0 


g  (x+A  x q )  +  (p(Tn-eA)k(Tn-6A,x)A 


'0 


0 


Hb 


•0 


g(x'4A-TQ)  +  cp(Tn-eAA)k(Tn-eAA,x«  )A 


+  /  9(y)k(y,x> )dy 


K(x,x»;Tn,T  )  - 


'0 


q>'1*  ~J  K(x,y;r0,T1)k(y,x')dy 

T 


dx1 


0 


+  d?(A) 


=  g(x-{A-T0)  +  9(T0-eA)k(T0-eA,x)A 

+  J1  K(x,xt;T0,T1)g(x,4A-T0)dxt 
T0 

+  v('V0~6^A  f  1  K(x*yj 

To 

+  {9  (A ) . 


The  above  result  coincides  with  Eq.  (3.6). 
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